Dependency of the drag coefficient on boundary layer stability beneath drifting sea ice in the central Arctic Ocean

The ice-ocean drag coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{w}$$\end{document}Cw and turning angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta_{w}$$\end{document}θw are crucial parameters in ice-ocean coupled simulations, determining the transfer of momentum between the two media. These parameters are often treated as constants regardless of the static stability at the ice-ocean interface. This study investigates the variability of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{w}$$\end{document}Cw and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta_{w}$$\end{document}θw based on direct observations of thermal and kinetic energy balance. The observations were conducted beneath multiyear ice packs widely across the central Arctic during a period transitioning from ablation to refreezing, indicating significant variability of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{w}$$\end{document}Cw = 1–130 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document}× 10−3 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta_{w}$$\end{document}θw =  − 19–1° at 5 m depth. Comparing different stations, the observations suggest a pronounced dependence of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{w}$$\end{document}Cw on the stability parameter (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}μ) resulting from mechanical and buoyant forcing. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{w}$$\end{document}Cw rapidly decays with increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}μ, indicating that the ice-to-ocean momentum transfer is enhanced for neutral or unstable conditions, while it is weakened for stable conditions. In addition, observed vertical profiles of currents revealed that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\theta_{w}|$$\end{document}|θw| tends to be smaller for unstable and larger for stable conditions. We suggest that numerical simulations using constant values could result in an underestimate of large-scale near-surface currents during the ice growing period.

www.nature.com/scientificreports/need to reevaluate those parameters in the context of the ongoing global warming and diminishing sea ice extent in the Arctic Ocean 17,18 .
For the bulk estimate of interfacial turbulent fluxes, particularly when using data from autonomous drifting buoys, the differential velocity between the ice and ocean is commonly neglected [24][25][26][27] .In such instances, the stress is frequently parameterized solely by the sea ice drift ( U i ) with some empirical constants, e.g., τ w = 0.0104U 1.78 i 24,25 .Alternatively, it can be determined by solving the implicit equation of Rossby similarity assuming a neutral IOBL 27 : U i u * 0 ∝ ln u * 0 fz 0 , where u * 0 = √ τ w is the interfacial friction velocity, f is the Coriolis frequency, and z 0 is the surface roughness.
In the Arctic Ocean, there are relatively few reports of the drag coefficient based on direct measurements of the Reynolds stress e.g., 11,12,[28][29][30][31][32] .For example, Shirasawa and Ingram 30 proposed a value of C w = 5 × 10 −3 at 1 m beneath the ice bottom in the Barrow Strait of the Canadian Arctic.Some numerical models solve the 'law of the wall' (LoW) problem for the determination of C w assuming a logarithmic current profile in the vertical direction under a neutrally stratified IOBL, i.e., , where κ is the von Kármán constant, and z is the vertical coordinate and points upward.For z = − 2.0 m, the above model yields an estimate of C w = 12 × 10 −327 .Cole et al. 12 performed year-long measurements of C w and θ w utilizing autonomous drifting buoys in the marginal ice zone of the Canada Basin.They found seasonal and regional variations in these parameters, particularly noting an anomalously reduced C w during the melting season compared to the rest of the year.Heorton et al. 11 evaluate those parameters in inverse modelling techniques based on the momentum balance of freely drifting sea ice, proposing an average value of C w = 2.4 × 10 −3 .
From the viewpoint of the steady-state momentum balance of the sea ice, the drag coefficients may significantly affect the ice drift, as described by the relation below 9,33 : where Na = √ ρ a C a /ρ w C w is the Nansen number; U a is the surface wind, C a is the air-ice drag coefficient; ρ a are ρ w are the densities of air and water, U wg is the geostrophic current, and θ a is the turning angle in air.We assume that the ice internals tress is ignorable compared to the other kinematic terms, resulting in the state of free drift 9 .In Eq. ( 2), Na can be regarded as the wind factor representing the ice drift relative to the surface wind speed [34][35][36] , assuming a negligible geostrophic current, i.e., U wg ≈ 0. If the stability state differs between the fluids below and above the ice layer, the ice drift speed could vary significantly through Na.
This study aims to derive the horizontal distribution of C w and θ w across roughly 2500 km in the central Arctic Ocean using measurements obtained during the RV Polarstern expedition PS138 in 2023 (Fig. 1); 37,38 .We explore the relationship of the key variables C w and θ w against the static boundary condition at the ice-ocean interface and address the question of how the static stability parameter related to basal melting and refreezing can affect the exchange of momentum across the ice-ocean interface and ultimately influence the large-scale circulation in the central Arctic basins.

Results
We conducted under-ice observations at nine ice stations during the PS138 cruise in the Eurasian Basin, encompassing the Amundsen and Nansen Basins (see Methods for further details; Fig. 1a).The observations were based on the eddy-covariance technique, which directly measured the three-dimensional movement of seawater in a small sampling volume beneath drifting ice floes.Based on these data, we calculated the cross-spectral power for the turbulent fluxes of kinematic and thermal energy, respectively represented by τ w = �u ′ w ′ � 2 + �v ′ w ′ � 2 and F h = ρ w C p �w ′ T ′ � .In the following sections, we will display the results of C w , which is calculated from the relation in Eq. (1), based on the ECS measurement of τ w and the ice-referenced velocity U 0 .Furthermore, the turning angle ( θ w ) is estimated from vertical profiles of horizontal current based on an ice-fixed acoustic Doppler current profiler (ADCP) (see the corresponding subsection in the Methods).

Turbulent fluxes from ECS observations
For Stations #1 to #5, all turbulence-related variables obtained by the eddy-covariance system (ECS) are well correlated with the mean current intensity U 0 (equivalent to U w − U i ) directly measured from the drifting sea ice platform (Fig. 2; Figs.S1-S9).These include turbulent kinetic energy , and turbulent heat flux ( F h ), where ( u ′ , v ′ , w ′ ) is the eddy component of the ocean current observed in the ice-fixed frame, T ′ is the eddy component of the seawater temperature, ρ w the seawater density, and C p the specific heat.For example, in the case of Station #3, the time series of burst observation confirms a clear concurrency between U 0 and the low-passed ice drift ( U i ), which were measured by independ- ent instrumentations (Fig. 2a,b).Furthermore, Q and u * 0 exhibit corresponding temporal variations with U 0 and consequently U i (Fig. 2b,c).In other words, the faster the ice floe drifted, the more energetic the turbulent mixing became.At Station #3, the variation in F h is well explained by U i , rather than by the change in local temperature immediately underneath the ice bottom (Figure not shown).
The friction velocity ( u * 0 ) varied moderately from station to station along with U 0 , as shown in Fig. 3b.Mean- while, F h exhibited pronounced changes in time and location.Results from the earlier stations, particularly Stations #1, 3, and 5, reveal that a substantial amount of heat was transferred from the ocean to the ice bottom through turbulent mixing.In particular, the results from Station #4 show a median estimate of F h = 15 W m −2 with an interquartile range of 30 W m −2 , when the floe was in a region with a discernibly reduced ice concentration (Fig. 1a).During those stations, the far-field water temperature ( T w ) was above the freezing point ( T f ), determined by the practical salinity ( S w ) at the same depth, where the gap between T w and T f fell within the range between 0.2 and 0.4 K (Fig. 3a). (2) In contrast to the earlier stations, Stations #6 to #9 exhibited extremely small or nearly zero F h (Fig. 3c).This can be explained by the presence of near-freezing seawater in the upper part of mixed layer ( T w ≈ T f ) (Fig. 3a).Especially for Station #9, the observed F h was mostly negative, where turbulent heat flux was directed downward at the interface, i.e., from the ice to the ocean.
The ECS observations exhibited that the medians of the drag coefficient C w largely increased with increasing station number, from the value of C w = 1 × 10 −3 to 130 × 10 −3 , transitioning from early August to late September (Fig. 3e).Meanwhile, the turning angle θ w generally diminished in terms of its magnitude, e.g., from θ w = − 19 • to 1 • at 5 m depth, with the exception of the relatively small deviation of θ w = − 5 • at Station #1 (Fig. 3f).www.nature.com/scientificreports/

Conductive and latent heat
The SIMBA ice mass balance buoys (refer to Methods) demonstrated the temporal variation of ice thickness ( h i ) as the ice floes drifted (Fig. 4b).For example, Station #1 exhibited h i decreasing by more than 0.4 m for three weeks in August until it underwent a transition to the growing phase around 10th in September (inverted triangles in Fig. 4c).Stations #4 and #8 switched to the growth in mid-October, while Station #5 did that in late October.At the remaining stations of #6 and #7, h i remained relatively constant through the end of the year.
According to the time evolution in h i , detected by tracking the SIMBA's ice bottom, the latent heat ( L f ) can be estimated (refer to Methods).For Stations #1 and #4, L f exceeded as large as more than 50 W m −2 in nega- tive values during the melting period (Fig. 4c).For Stations #6 and #7, deployed in the highest latitudes, L f is estimated to be nearly zero through the year.The rest of stations (#5 and #8) showed reduced variability in L f being less than 10 W m −2 in magnitude.
At Stations #2, 3 and 9, there are no records of T i available for the calculation in terms of q .We hence estimated q by assuming a linear vertical profile of T i , between the measured near-surface air temperature ( T a ) and the freezing-cold temperature in water, i.e., T w ∼ T f (Fig. 4a) (see Methods).Using this approximation, we derived the rough estimates of q = − 2.0, − 1.5 and 4.4 W m −2 , respectively for Stations #2, 3 and 9 (Table 1).In panels (b-d), red lines show mean current intensity ( U 0 ) directly measured by ECS, while colored lines represent ADCP measurements at levels of 3, 7 and 12 m (see legend for corresponding depths).In (d), magenta vertical bars denote negative F h .
For Stations #1 and #4, the measured F h values were so small (10 and 17 W m −2 , respectively) that they can- not explain the rapid decay of the ice column at the undersurface (Table 1).The underrated F h based on the ECS would be attributed to some environmental factors.Oceanic heat content stored in the fresh meltwater, which lies between the ice and the ECS instruments, has the potential to facilitate basal ablation 39 .

Static stability near the interface
We now address the static stability issue drawing from the measured heat and momentum balance.To characterize the stability state in IOBL, the introduction of the Obukhov length scale as is advantageous, where �w ′ b ′ � represents the turbulent buoyancy flux at the interface.The buoyancy flux ( �w ′ b ′ � ) is estimated based on thermohaline conservation at the interface, utilizing the hydrographic observations for the far-field temperature and salinity (refer to Methods) as well as the ice mass balance observations for the conductive heat flux.In (a), a black dashed curve denotes an average of T a from the six SIMBA buoys, while in (d) it denotes q derived from the averaged T a and the constant ice thickness of 1.2 m.In (c), inverted triangles denote zero-crossing points in L f .
In general, �w ′ b ′ � concurrently varied with F h (Table 1).For the early stations (Stations #1 to #5), which occurred in August, the ice-ocean interface exhibited �w ′ b ′ � > 0 , likely due to the significant melting of sea ice at its bottom, resulting in the stabilized boundary condition.Conversely, for the later stations (#6 to #9), �w ′ b ′ � was estimated to be significantly lowered.Especially, for Stations #7 and #9, �w ′ b ′ � indicated slightly negative values.Thus, the interfacial boundary condition can be perceived as transitioning from stable to neutral and/ or even slightly unstable.
We then transform Lo into a non-dimensional form, known as the stability parameter ( µ) 27 (Fig. 3d): where µ represents the ratio of the planetary boundary layer thickness 4) is formu- lated to draw an analogy with the physics of the air-ice boundary layer 40,41 .General interpretation of physics can describe that µ > 0 indicates a stable boundary condition, typically associated with ice ablation, while µ < 0 indicates an unstable condition with brine rejection due to ice growth.Hence, µ = 0 signifies a neutral boundary flux condition.Using the observed u * 0 (Table 1), we also note that the mean value of L PB can vary between 110 and 580 m.According to the observations, the earlier group (#1-5) in August typically exhibited positive values of µ ranging between the orders of 0.1 and 1, indicating a stable boundary condition.In the meantime, the latter group (#6-9) in September showed nearly zero or slightly negative values of µ , respectively indicating neutral or statically unstable boundary conditions (Fig. 3d).

Vertical structures of Ekman currents
The vertical structure of horizontal currents in IOBL was examined using data from an ADCP (Fig. 5).The measured current speed ( U 0 ) relative to the drifting sea ice exhibited a characteristic vertical variation of the LoW theory, logarithmically increasing with depth 27 : The current strength increased with depth and reached an equilibrium level around a depth of 10 m, where the geostrophic current ( U geo ) was expectedly dominant (Fig. 5a).For the normalized vertical profile of Û0 = (U 0 (z) − U geo )/U geo , where U geo is derived as an average of U 0 (z) over a depth of 15 to 20 m, the surface roughness of |z 0 |= 4 cm 27 for the logarithmic model above is in good agreement with the observed vertical profile (Fig. 5b).
Comparing current speeds at different levels as detected by the ADCP (colored curves in Fig. 2b-d), it is evident that the ECS instrument was positioned in the middle of the LoW boundary layer, where the ice-relative current strength logarithmically increases along the depth (Fig. 5b).The current observed at the ECS level, typically 0.7 m below the ice bottom, was as small as 0.3 to 0.6 times the current at 12 m depth.
In the present study, the ADCP observations reveal that the horizontal current vector was rotated clockwise relative to the surface current and mostly increased linearly with depth, as expected from the classical Ekman theory 5 (Fig. 5c).From Stations #1 to #8, the observation show rotation angles ranging from θ w = − 18.8° to − 3.5° and from θ w = − 26° to − 6.4°, respectively at 5 and 10 m depth (Fig. 3f).We note that θ w is a deviation from

Discussions Dependency of C w and θ w on stability parameters
Based on Eq. (1), the relationship between the ice-referenced current ( U 0 ) and the Reynolds stress ( τ w ) is investi- gated for all ice stations (Fig. 6).For every ice floe, there was a clearly positive relation between the two variables, where the stress τ w displays a quadratic growth with the mean relative current U 0 (Fig. 6a-i).
The changing rate of C w as a function of U 0 exhibits considerable variability among different ice stations (Fig. 3e).For Stations #1 to #5 in mid-to late-August, the data points are distributed to fit the curves of C w = 1 to 5 × 10 −3 (Fig. 6a-e).These values for C w appear slightly smaller yet but remain within the same order of magnitude as those derived from some previous studies employing direct observations using the eddy covariance method (c.f., C w = 5.4 × 10 −3 by Shirasawa and Ingram, 1997 30 , horizontal line in Fig. 7).In contrast, for Stations #6 to #9, τ w exhibits a higher rate of increase with U 0 , indicating a significantly larger value of C w ranging approximately from 10 to 100 × 10 −3 (Fig. 6f-i).
By considering the median values, the universal characteristics of C w derived from the ECS are extracted in the context of their dependence on the turbulent heat flux ( F h ) and the non-dimensional stability parameter µ (Fig. 7).The scatter diagram exhibits a pronounced feature of the C w curve, notably decaying with increasing values of F h and µ.
For the analogy of the ice-air boundary layer physics 41 , C w can be expressed as a function of µ (blue curve in Fig. 7b): where z is taken as the planetary boundary layer thickness, i.e., |z| = L PB and 340 m; the roughness length |z 0 | is 4 cm.The similarity functions of A(µ) and B(µ) are defined for the stable ( µ ≥ 0 ) and the unstable domains ( µ ≤ 0 ), respectively by the linear and nonlinear (the reciprocal of cube root) functions in terms of µ 41,42 : and This interpretation suggests that for the stable condition, i.e., µ ≥ 0 , the drag coefficient C w should decrease relatively slowly with increasing µ (Fig. 7b).Regardless, for the unstable condition of µ ≤ 0 , the value of C w significantly grows as µ decreases across the neutral point of µ = 0. Our results also suggest that the melting condition strengthens local stratification near the interface, thereby hindering the momentum transfer due to reduced drag efficiency.On the other hand, a transition toward basal ablation disrupts near-interface stratification, promoting vertical momentum transfer due to increased drag efficiency 24 .
Based on their direct measurements of turbulent momentum fluxes ( τ w ) at approximately 2.5 and 4.5 m below the ice bottom, Cole et al. 12 reported that C w varied in an approximate range of 1-20 × 10 −3 for the weekly median values, with the minimum occurring during the melt season.Qualitatively, the temporal variation suggested by ( 8) www.nature.com/scientificreports/this study is consistent with their results.However, for more quantitative discussions, our estimate ( C w > > 10 × 10 −3 at µ ≤ 0 ) is quite large, perhaps due to the relatively short observational period, typically lasting for a few days for each station.
According to the ADCP observations, for Stations #1 to #8, the turning angle ( θ w ) exhibited the certain rela- tionship against the parameters of F h and µ .Unfortunately, the statistical dispersion of an interquartile range for θ w (indicated by error bars) is too wide, typically exceeding 10°, to argue for a universal curve for θ w as a func- tion of F h and µ (Fig. 3f; Fig. 8).However, it is evident that the clockwise rotation by θ w tends to increase with increasing F h and µ for more stable conditions (Fig. 8).This can physically be interpreted that the Ekman flow rotates more quickly with depth in stable stratification, whereas it rotates more slowly in unstable stratification 10 .
The ice-air boundary layer theory 41 proposes the solution of turning angle θ w as a function of µ , based on the similarity functions of A(µ) and B(µ) 42 (blue curve in Fig. 8b): The above semi-empirical model largely aligned with the observed variation of θ w along the coordinate of µ .The solution also provides diagnostic solutions of |A|= 11.0 and B = 0 for the neutral point ( µ = 0).These show a certain discrepancy from the widely-used approximations of |A|= 2.3 and |B|= 2.1 27 , presumably due to the distor- tion caused by the extreme values at µ ≤ 0 observed in our data.The validity of the presented constants will be verified after accumulating knowledge from unstable-time surveys.

Analogy with atmospheric boundary layer
The dependency of the oceanic drag coefficient and turning angles on the static stability presents a strong similarity with the ones in the atmospheric boundary layer in pack ice 9,40,41 .In the stratified air boundary layer, C a tends to be extremely small, typically O(10 −4 ), and the turning angle from the geostrophic wind is typically as large as θ a = 15° counterclockwise in the Northern Hemisphere.Meanwhile, in unstable conditions like in the marginal ice zone 43 , C a increases to be 10 times larger, while the rotation decreases approaching θ a = 0°.

Influences on sea ice movement
In the multiyear pack ice during mid-summer, warm water in melt ponds or open water in leads is heated by solar radiation.This process potentially weakens the stratification within the air boundary layer, leading to statically unstable conditions.Conversely, during the same period, the oceanic boundary layer remains stable due to the prevalence of fresh meltwater near the surface.In the disparate situation, the value of Na = √ ρ a C a /ρ w C w (see Eq. 2) could grow strongly due to the combination of increasing C a and decreasing C w .Indeed, it is acknowledged that the wind factor ( α ) varies seasonally by 50% or more between the summer months (July to October) and the winter months (December to March) in the central Arctic 36 .The seasonal variation in α is often attributed to reduced (increased) internal stress among ice floes during summer (winter) 34,35 .Contrary to the general belief, the findings from the present study suggest that the disparate stability states between the air and oceanic boundaries could partly explain the seasonal changes in α ≈ Na.

Implication to ice-ocean coupled simulations
Changes in the parameterization of ice-ocean drag play a significant role in ocean dynamics, affecting processes such as seasonality and trends in ocean surface stress, and the dynamics of the mixed layer 17,18 .Sensitivity studies have demonstrated that different parameterizations introduce significant spatial and temporal differences in the distribution of these quantities over the Arctic Ocean compared to models using constant drag coefficients 14 .Most studies replacing the commonly used constant drag coefficient with more accurate representations of iceocean friction have shown improved model results 14,17,18 .Primarily, these papers discuss form factors and various parameterizations of friction between the ocean and ice.
From the perspective noted above, the findings from this study provide an important instance of the potentially variant drag effect from the real field.Our results suggest that numerical simulations with invariant C w could significantly underestimate the kinetic energy exchange between the drifting pack ice and the underlying ocean, especially when the ice is growing.This would lead to potential misrepresentation in the simulated largescale circulation, the sea ice drift and the near-surface currents.We expect that the incorporation of stabilitydependent C w could lead to an improvement in the predicted volume export of sea ice through the Fram Strait, which may partly control the pan-Arctic sea ice extent 4,7 .www.nature.com/scientificreports/Technically, the ice-ocean coupled model can incorporate the empirically derived relation (Eq.6) to update those values based on predictions of the stability parameter µ (Eq.4) from the interfacial buoyancy flux w ′ b ′ and friction velocity ( u * 0 ).An estimate of w ′ b ′ requires solving the heat and salt balance around the ice-ocean interface.The calculation should be based on variables averaged over a period of 1 day or perhaps longer to avoid frequent updates in the calculation.
We also emphasize that the relationship for the stability-dependent C w and θ w were derived from a limited dataset, and therefore, we recommend that the coefficients of the similarity functions (Eq.7-9) be updated and evaluated in terms of the simulated ice-drift properties and the eventual sea ice distribution.

Concluding remarks
In this study, we conducted the direct observations of turbulent fluxes of momentum and heat near the ice-ocean interface as well as vertical profiles of ocean currents and sea ice mass balance in the central part of the Arctic Ocean.The observations were conducted at nine ice stations across the Amundsen and Nansen basins as part of RV Polarstern expedition PS138 (ArcWatch 1) in August/September 2023 (Fig. 1).
The eddy-covariance measurements refined the estimate of the drag coefficient ( C w ) at the interface, reveal- ing a greater variability than previously assumed, ranging within C w = 1 to 130 × 10 −3 (Fig. 3).Our observations demonstrate a strong dependence of C w on the stability parameter µ = L PB L O (Figs. 6, 7).Specifically, C w exhibits an exponential decay as the flow becomes more stably stratified by basal ice melting, and conversely it increases as flows become neutral or slightly unstable due to ice growth.Additionally, from the ADCP observations, the turning angle ( θ w ) indicated such spatiotemporal variability as a function of µ , with medians ranging from θ w = − 18.8° to 1.3° at 5 m depth for the range of µ = − 2 to 5 (Figs.5-8).Overall, the degree of clockwise rota- tion ( θ w < 0 ) tends to increase for the stable condition, while it decreases for the neutral or unstable condition.
In the present study, the in-situ data from unstable conditions were limited to the early transitional period during the late summer in the central Arctic, when the cold front of the ice temperature did not reach the iceocean interface yet (Fig. S10c-h).Future observations should aim to cover a wider range of variability for C w and θ w regarding the stability parameter, especially focusing on the negative domain of µ during the mid-winter.

Study domain
All observations presented here were obtained as part of the research expedition 'PS138/ArcWatch-1' 37,38 to the Arctic Ocean onboard the German ice breaker RV Polarstern (Alfred-Wegener-Institut; 44 .During this expedition, nine ice stations (Stations #1 to #9) were conducted in the Nansen and Amundsen Basins between 8 August and 19 September 2023 (Fig. 1a).
Station #1 was located at (84.07°N, 31.21°N)northeast of Svalbard, and Stations #2 to #4 were located farther east at the similar latitudes (see Table 1 and Fig. 1a for precise locations).Stations #5-#7 were located on ice floes within the approximate meridional range of 110°-130°E, where Station #7 was positioned near the geographical North Pole.Stations #8 and #9 were positioned along the 60°E line, respectively at latitudes of 88.0°N and 85.5°N.

Ice thickness variation
Ice thickness observations based on helicopter-borne electromagnetic (EM) induction sounding 45 were repeatedly conducted during the expedition.The airborne surveys show that the mean ice thickness statistically fell within a range of 1.3-1.4m for the observational areas (personal communication, J. Rohde).Snow thickness atop the ice was typically around 0.1 m or less.
The ECS instrument was selectively deployed on undeformed part of each ice floe, where the thickness ranged approximately between 1.2 and 1.5 m.It is noted that at Station #1, the ice floe was uniformly thin, measuring approximately 0.75 m throughout.

Turbulent flux measurements
Measurements of momentum and heat turbulent fluxes were collected using an eddy-covariance technique [46][47][48] during all nine ice stations (Fig. 2).The ECS consisted of a 600 MHz Nortek Vector (herein referred to as Vector) three-dimensional current meter (Nortek, Norway) combined with a highly accurate RINKO-EC ARO-EC-CM (herein referred to as RINKO-EC) temperature probe (JFE Advantech Co., Japan).The ECS was deployed through a hydrohole in the ice, and the Vector transducer was positioned approximately 0.5 to 0.8 m beneath the ice bottom.The spearhead of the RINKO-EC sensor was positioned within 10 mm from the lateral boundary of the sampling volume (diameter = 14 mm; height = 14 mm), located about 157 mm from the central transmitter, and inserted at an angle of 45 degrees from below relative to the vertical axis of the water current measurement by the Vector.
The Vector sampling rate was set to 8 Hz during 256-s burst intervals.The ECS collected 2048 data samples per burst, and a burst was set to be recorded every 10 min.The current velocities ( u, v, w ) were measured in the instrument-fixed Cartesian coordinate ( x, y, z ), where a positive z denotes an upward direction.The Vector cur- rent velocity is stated to have a nominal precision of ± 1 mm s −1 or ± 0.5% for the measured value.The RINKO-EC response time for the temperature is stated to be < 0.3 s 48 .The sensor was calibrated at the manufacturer prior to the expedition, ensuring a precision of ± 0.002 °C.
Prior to the spectral calculations, the raw data underwent conventional postprocessing.Attitude information recorded onboard such as heading, pitch and roll were utilized for the correction of inclination within 5 degrees.Data outside these criteria were discarded.Outliers beyond three standard deviations were removed, followed by interpolation with a linear function 31 .Subsequently, the deviations were extracted by demeaning the sequential records.The cross-power spectral density (CPSD) was calculated using P. Welch's periodogram method 49 , where where ω represents the frequency, and S wC denotes the CPSD between vertical flow w and an arbitrary variable C.

Vertical profiles of currents
For all ice stations, horizontal ocean current velocities were obtained using a Nortek 1 MHz Signature1000 ADCP 51 (Fig. 5).The ADCP was deployed through a hydrohole within 5 m from the ECS location as the transducer looks downward.Three-dimensional current velocities were collected at a sampling frequency of 10 min and a burst interval of 1 min, with a sampling rate of 1 s during each burst.The vertical extent of the current vector typically ranges from 1 to 20 m from the ice bottom, with a vertical bin size of 1 m and total of twenty bins.Current velocity was recorded in the instrument-fixed Cartesian coordinate system (XYZ).During postprocessing, horizontal velocities with a Percent Good ≤ 60% were discarded.The reported precision of the ADCP is ± 0.3% of the measured current speed.

Far-field temperature and salinity
Far-field temperature ( T w ) and practical salinity ( S w ) are calculated as averages of respective hydrographic param- eters over 5-10 m depth using profiles obtained by a tethered, free-falling MSS-90L microstructure profiler (Sea and Sun Technology, Germany; 52 (red filled circles in Fig. 1b).Typically, 20-40 profiles were obtained during each ice station.The MSS was lowered through a hydrohole in the ice, located ~ 30 m from the ECS, and descending to a maximum depth of 300-400 m.Sensors mounted on the instrument included seawater conductivity (SST small), temperature (PT100) and pressure (PA7-100) (CTD) with a precision of ± 0.002 mS/cm, ± 0.002 °C, and ± 0.1 dbar, respectively.The response time of these sensors is approximately 0.15 s, resulting in a vertical resolution of 0.05-0.10m for the raw CTD variables.The MSS data are preliminary, with an accuracy of 0.05 in practical salinity determined by comparison to data from more accurate shipboard CTD (not shown).This is sufficient for the analysis in this work.The raw data was processed using the MSP toolbox developed by the Institute of Baltic Research (Warnemünde, Germany; 53 .Further observational protocols and data postprocessing related to the MSS-CTD data are described in more detail in Kawaguchi et al. 54 .

Ice mass balance observations
Ice interior temperatures ( T i (z) ) were measured by SIMBA-type ice mass balance buoys equipped with ~ 5 m long thermistor strings (SAMS Enterprise, Scotland; 55 (Fig. 4; 56,57 .The SIMBAs were deployed in undeformed ice on Stations #1, 4, 5, 6, 7, and 8, at a distance within 100 m from the ECS, where initial thickness of the ice was respectively 0.8, 2.0, 1.3, 1.5, 2.7, and 1.3 m, respectively (Fig. 4b).The SIMBAs measured T i at a vertical resolu- tion of 0.02 m, covering the depth range from the ice surface to 4.8 m.The sampling interval was set to every six hours, and the manufacturer reports an accuracy calibrated to ± 0.125 K.The ice-ocean interface is visually identified from a distinct gap between the two media based on 30-s and 120-s heating mode temperatures of SIMBAs (green lines in Fig. S10; 58 .

Thermohaline balance at the interface
One of the main aims of this study is to investigate the ice-ocean drag coefficient in terms of the stability in IOBL (blank red circles in Fig. 1b).However, the ECS had no sensor to measure salinity flux, and consequently it did not allow a direct estimate of the buoyancy flux �w ′ b ′ � or any parameter related to the interfacial stability.Instead, we tried to deduce �w ′ b ′ � based on the estimate of turbulent temperature flux �w ′ T ′ � from the ECS observation.
For that purpose, we first used the following approximation for the heat balance 57,58 : where L f (t) is the latent heat, q(t) is the conductive heat flux, t is time, h i is ice thickness, ρ i and ρ w are the den- sity of ice and seawater and 910 and 1023 kg m −3 , respectively; the constant for latent heat ( Q L ) is 276 kJ kg −1 .Using the autonomous ice mass balance observations by the SIMBAs 56 , q(t) is calculated from the vertical ice temperature gradient near its bottom 59,60 according to where T i (z) is the SIMBA-based ice temperatures along the z-coordinate.We calculate an average value of dT i (z) dz over a vertical extent of 0.2 m (equivalent to 10 thermistors) above the ice bottom.The thermal conductivity of sea ice ( K i ) is approximated by 61 : where K f is the thermal conductivity for fresh ice (2.04 J m −1 K −1 s −1 ) and S i is the salinity of sea ice and assumed to be 6.0.In the equation above, the ice temperature ( T i _btm ) near the bottom is assumed to be the freezing point of seawater ( T i _btm = T f ≈ − 1.8 °C).From the SIMBA observations, q is typically directed downward and estimated to be − 6.0, − 1.0, − 1.8, − 4.2, 0 and − 4.0 W m −2 respectively for Stations #1, 4, 5, 6, 7, and 8 (Table 1); then it is substituted into Eq.( 12).At Stations #2, #3 and #9, there are no ice temperature profiles available.Thus, we estimated q according to a linear vertical profile in T i (z) deducing the temperatures near top and bottom of the ice column.Regarding the near-bottom temperature, it was set ad hoc to the freezing point of T f ≈ − 1.8 °C for average seawater [57).At those stations, the air temperatures were deduced from the mean timeseries profile constructed from the remaining SIMBA temperature data (black dashed curve in Fig. 4a).In the calculation, the ice thickness was assumed to be constant at 1.2 m.
The salinity flux �w ′ S ′ � at the interface is derived from the salt balance, in which �w ′ S ′ � is proportional to the melt rate at the ice bottom 27 : where S 0 is salinity at the interface, and it is derived by solving the following quadratic equation 27,54 : where, α h and α s are the extraction coefficients for temperature and salinity, respectively; m is the slope of 'freez- ing line' for the freezing temperature from the UNESCO formula 62 as a function of salinity: In the calculation, the certain combination of α h α s = 50 and α h = 0.0055 are set considering the double-diffusive effects from fresh meltwater accumulated near the ice-ocean interface 27,31,63 .
The buoyancy flux �w ′ b ′ � is obtained by an equation 27 below: where β S are the haline contraction and β T thermal expansion factors evaluated at T 0 and S 0 following the method of McDougall 64 .In Eq. ( 19), we directly substituted the readings of �w ′ T ′ � from the ECS.

Figure 1 .
Figure 1.(a) Map of the study area, including sea ice concentration on September 15, 2023.Red dots represent Stations #1-9.Gray contours show the IBCAO bottom relief:100, 500, 1000 and 2000 m.(b) Potential temperature and salinity of sea water at each station, at depths within mixed layer (filled circles; T w and S w , respectively) and the ice-ocean interface (blank circles; T 0 and S 0 , respectively).Contours show the potential density anomaly.

Figure 2 .
Figure 2. Time series of turbulent fluxes from the ECS observations at Station #3: (a) ice drift intensity (thin line: unfiltered; bold line: 12-h lowpass filter); (b) turbulent kinetic energy ( Q ); (c) friction velocity ( u * 0 ); (d) turbulent heat flux ( F h ).In panels (b-d), red lines show mean current intensity ( U 0 ) directly measured by ECS, while colored lines represent ADCP measurements at levels of 3, 7 and 12 m (see legend for corresponding depths).In (d), magenta vertical bars denote negative F h .

Figure 3 .
Figure 3. Overview of statistics for kinematic and thermodynamic parameters during all stations: (a) temperature deviation from the freezing point ( T = T w − T f ), (b) friction velocity ( u * 0 ), (c) turbulent heat flux ( F h = ρC p �w ′ T ′ � ) in black and conductive heat flux ( q ) in gray, (d) stability parameter ( µ = L PB /L O ), (e) drag coefficient ( C w ) and (f) turning angle ( θ w ) at 5 m depth.Bold black bars indicate medians with red error bars representing first-to-third quartiles.Numbers in red circles at the top denote the corresponding station.In (b), blue line shows mean current intensity, U 0 , divided by 10.

( 3 )Figure 4 .
Figure 4. Time series of SIMBA-derived variables: (a) air temperature T a , (b) ice thickness h i (c) latent heat L f and (d) conductive heat flux q .In (a), a black dashed curve denotes an average of T a from the six SIMBA buoys, while in (d) it denotes q derived from the averaged T a and the constant ice thickness of 1.2 m.In (c), inverted triangles denote zero-crossing points in L f .

Figure 5 .
Figure 5. Vertical profiles of medians of horizontal current velocity from Stations #1 to #9: (a) current magnitude U 0 (z) , (b) its normalized deviation ( Û0 (z) = U 0 − U geo /U geo ), where U geo is the geostrophic current averaged over depths of 15-20 m, and (c) turning angle θ w relative to the near-surface current, where positive values indicate counterclockwise rotation, and (d) number of quality data where Percent Good ≥ 60%.In (b), a black dashed curve represents the LoW logarithmic model, with a surface roughness of z 0 = 4 cm.

Figure 6 .
Figure 6.Relationship between the current speed ( U 0 ) measured from the ice-fixed frame and the ice-ocean Reynolds stress ( ρ w τ w ) for each ice station.Black curves show the drag coefficients ( C w × 10 3 ) for the quadratic model (Eq.1).

Figure 7 . 2 0.
Figure 7. Dependence of C w on (a) F h and (b) µ , where µ = L PB L o = κ�w ′ b ′ � fu * 2 0 .Black dots indicate medians for each station, with first-to-third quartiles shown by horizontal and vertical error bars.Bold blue curve shows the least-square regression deduced from the observed data.Note that the vertical axis is shown on a log scale.Gray horizontal lines represent the constant suggested by Shirasawa and Ingram 30 .

Figure 8 .
Figure 8. Same as Fig. 7 but for the turning angle θ w .Bold circles and triangles denote the values at 5 and 10 m from the ice bottom.In (b), a blue curve shows the solution of θ w (µ) between the near-surface current and the geostrophic current, constructed with the regression coefficients for C w (µ).

Table 1 .
Overview of under-ice observations during the Polarstern ArcWatch expedition, 2023.Turning angles ( θ w ) are shown for depths of 5 m and 10 m in the top and bottom rows, respectively.Numbers inside and outside parentheses denote medians and interquartile ranges, respectively.